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SECTION  I 


INTRODUCTION 

Time-dependent  computations  have  been  used  primarily  as  a 

device  to  obtain  steady  state  solutions  of  gas  dynamical  problems. 

As  such,  they  do  not  differ  conceptually  from  the  relaxation 

1,2 

method  ;  no  computed  results,  except  the  final  ones,  are  to 
be  interpreted  physically.  In  principle,  one  can  use  a  numerical 
technique  which  is  consistent  with  the  physico-mathematical  model 
only  in  the  limit,  when  all  time  derivatives  vanish  (see,  for 
example.  Refs.  3  and  13). 

In  the  first  impressive  applications  of  the  relaxation 

method  #  the  computations  (performed  by  hand)  were  constantly 
kept  under  the  control  of  the  analyst  who,  perhaps  semi- :onsciously , 
would  use  his  physical  intuition  to  direct  the  wcrk  to  a  reasonable 
end.  lr  the  task  is  shifted  to  a  machine  which  operates, 
uncontrolled,  by  brute  force,  most  likely  the  results  do  not 
converge  or  converge  to  a  physically  unreasonable  pattern.  This 
can  be  said  of  relaxation  methods  and  time-dependent  techniques 
as  well.  So  long  as  no  physical  meaning  can  be  given  to  the 
time-dependent  results  as  functions  of  time,  physical  intuition 
cannot  be  used  to  explain  the  origin  and  nature  of  troubles.  Such 
vague  terms  as  "nonlinear  instability"  are  then  used  to  mean  that 
there  is  trouble  but  one  does  not  know  why. 

I  intend  to  shew,  by  a  detailed  analysis  of  a  particular 
problem  .that  time-dependent  techniques  may  be  devised  which 
closely  describe  an  actual  time-dependent  evolution;  that  errors. 
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generated  by  local  inconsistencies  of  the  numerical  model  with 
the  physical  model,  propagate  throughout  the  flow  field  according 
to  physical  laws;  consequently,  that  troubles  can  be  traced  back 
to  their  point  of  origin;  and  finally,  that  the  initial  and  boundary 
conditions  are  as  important  in  the  numerical  treatment  as  they 
are  in  the  physical  phenomenon. 

The  problem  chosen  for  analysis  is  the  two-dimensional, 
inviscid  compressible  flow  past  a  circular  cylinder.  The  choice 
ij  justified  by  some  interesting  features  of  the  flow,  namely, 

1)  The  cylinder  may  move  at  a  subsonic  speed.  In  this  case, 
the  steady  state  solution  should  cover  the  entire  plane. 

2)  Under  the  above  assumption,  a  local  bubble  of  supersonic 
flow  (relative  to  the  cylinder)  may  form.  In  this  case,  one 
expects  an  imbedded  shock. 

3)  If  the  cylinder  moves  at  a  supersonic  speed,  the  perturbed 
region  in  a  steady  state  is  limited  to  a  portion  of  the  plane 
(behind  ti  e  bow  shock)  . 

4)  Under  the  above  assumption,  can  a  numerical  solution  be 
found  which  has  all  the  characters  of  an  inviscid  flow,  that  is, 
without  separation  in  the  leeward  side? 

5)  How  does  a  steady  pattern  form  when  the  cylinder  takes  on 
a  constant  speed  starting  from  a  state  of  rest? 

In  addition,  the  choice  of  a  circular  cylinder  instead  of 
another  geometrical  shape  is  justified  by  some  simplification  in 
the  treatment  of  the  wall  points  (which  does  not  restrict  the 
general  conclusions) .  Finally,  a  similar  problem  has  been  con¬ 
sidered  by  other  authors^' and  myself^.  Thus. a  critical 
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comparison  can  be  made  which  confirms  some  of  the  conclusions  of 
the  present  paper. 

SECTION  II 

THE  PHYSICAL  PROBLEM 

Consider  a  cylinder  at  rest  in  a  gas  at  reot.  Let  t.he 
center  of  the  cylinder  star  moving  along  a  straight  line  (the 
x-axis),  from  right  to  left,  with  a  velocity  increasing  from  0 
to  a  maximum  value,  VQ.  For  example,  let 

(1)  x  =  -  V  sinuut  (0  <  t  <r  n/2uu) 

o  o  —  — 

where  t  is  time  and  x  is  the  abscissa  of  the  center  of  the 

o 

cylinder  on  the  fixed  x-axis.  For  t  >  tt/2uj  ,  let  the  cylinder 

move  at  a  constant  speed  along  the  x-axis.  Asymptotically  in 

time,  a  steady  flow  should  be  observable  in  a  frame  moving  with 

the  cylinder.  At  infinity,  such  flow  should  have  a  uniform 

velocity,  Vq,  and  uniform  values  of  the  other  physical  parameters 

(all  equal  to  their  values  in  the  gas  at  rest',,  We  want  to 

analyze  the  transient  from  a  physical  viewpoint  first. 

For  consistency  with  the  convention  adopted  in  the  numerical 

analysis,  nondimen sional  symbols  will  be  used.  All  pressures  and 

densities  are  scaled  to  their  values  (p _ ,  p_)  in  the  gas  at  rest, 

oo  oo 

all  velocities  are  seeded  ;o  Jd  To  #  all  lengths  are  scaled  to 

00  CD 

the  radius  of  the  cylindor,  r^,  and  all  times  are  scaled  to 
ro/  v'Pgp / .  In  particular,  the  speed  of  sound  in  the  gas  at 
rest  is  equal  to  in  nondimens ional  form,  where  y  is  the  ratio 

of  specific  heats.  A  pressure  coefficient  will  be  defined  as: 
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Our  arguments  will  be  illustrated  by  two  sets  of  sketches.  One 
will  depict  the  state  of  the  flow  at  a  given  time  on  the  plane  of 
motion,  by  an  isobar  pattern.  The  crotfs-section  of  the  cylinder 
will  appear  as  in  the  upper  part  of  Fig.  1.  The  other  will  show 
isobars,  as  functions  of  time,  along  the  heavy  line  in  the  upper 
part  of  Fig.  1,  which  consists  of 

1)  the  x-axis  in  front  of  the  cylinder, 

2)  the  upper  surface  of  uhe  cylinder,  and 

3)  the  x-axis  behind  the  cylinder. 

As  shown  in  the  lower  part  of  Fig.  1,  the  time  axis  runs  horizontally 
from  left  to  right.  Lines  AA  and  BB  represent  points  A  and  B, 
respectively,  in  a  frame  moving  with  the  cylinder. 

Between  BB  and  AA,  c=  6°/100- .9.  'f'his  strip  represents  the 
upper  part  of  the  cylinder.  Above  AA,  where  the  x-axis  in  front 
of  the  cylinder  is  represented,  £=-x- 0.1.  Below  BB,  where  the  x-axis 
behind  the  cylinder  is  represented,  Q-  -x+0.1.  In  a  steady  state 
(with  respect  tc  the  cylinder)  the  isobars  are  straight,  horizontal 
lines.  Any  perturbation  produced  by  the  cylinder  at  any  time  t  >  o 
is  confined  to  the  region  bounded  by  the  two  characteristics 
issuing  from  A  and  B  at  t=0,  unless  they  are  overcome  by  a  faster 
moving  shock.  The  upper  characteristic,  AM,  is  defined  by 
d£/dt=  yf“+  xq,  the  lower  one,  BN,  by  d£/dt=  -Jy"  +  xQ.  Therefore, 
once  the  cylinder  moves  at  a  constant  speed  V(=MQ,yY~)  both 
characteristics  become  straight  lines, defined  by 

o)  dc/dt-  -  yr~ [Mq  +  1), 

respectively. 
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Fig.  1 

When  the  cylinder  starts  moving  to  the  left,  the  portion  of 
the  cylinder  near  A  acts  on  the  flow  at  rest  as  a  piston  and,  in 
the  first  phase  of  the  motion,  the  ensuing  flow  is,  to  all 
practical  effects,  one-dimensional.  At  a  time, 

(4)  V  (y+1)$o<0)  |  =  (v+lfvj^ 

a  shock  builds  up  along  the  AM  characteristic.  Shortly  there¬ 
after,  the  one-dimensional  shock  reaches  its  peak  strength.  At 
tj-n/2uj,  when  the  cylinder  starts  moving  at  a  constant  speed,  the 
pressure  at  A  is  approximately  given  by: 


5 


(5) 


Jl. 

(1+  V  )  Y_1  M  1+  JT 
2  0 


if  V  «  1.  Thus 
o 


(6)  Cp^ti)  *2^Vo 

The  pressure  behind  the  fully  developed  shock  in  the  one- 
dimensional  case  is 


(7) 


P= 


2YM  »-(  Y-l) 

s 

Y+l 


where  M  is  the  ratio  between  the  modulus  of  the  shock  velocity 
s 

in  a  fixed  frame  and  the  speed  of  sound  in  the  gas  at  rest,  */yT 
The  velocity  of  the  gas  relative  to  the  shock  is 


(8)  urel=  v  +,y7  Mg 

where  V  is  the  absolute  velocity  of  the  gas  in  a  fixed  frame. 
In  front  of  the  shock  thus, 

urel=  YTM, 

Behind  the  shock. 


(9) 


urel=  YT  M 


(v-l)Msi  +  2 

s  ( Y+l )Ma 
& 


*  -  Ms 


Therefore, 


and 


2M  2  -  V  M  -2=0 

S  Jy  O  S 


(io>  v  i  { 3=  vo  +  v + « } 

By  substituting  M  from  (ID)  into  (7),  simplifying  for  V  «  1  and 

s  o 

comparing  with  (5'  or  (6)  ,  one  sees  that  the  final  pressure  between 
piston  and  shock  when  a  steady  state  is  reached  is  practically  the 
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same  as  the  pressure  at  the  beginning  of  the  constant  velocity 
phase  of  the  piston  motion,  if  VQ  is  small. 

At  B,  the  opposite  phenomenon  occurs.  If  B  were  a  one- 
dimensional  piston,  its  motion  would  produce  an  expansion  wave  and 
the  pressure  at  B  would  be  given  by 

lx. 

(11>  PB*  (1-  ■2^VV"1  *  1-/7'vo 

Consequently,  at  B 

us)  cp,r-2vr/vo 

We  can, thus,  anticipate  the  existence  of  a  first  phase  of  the 
motion  in  which  the  front  and  the  rear  of  the  cylinder  act 
independently  as  compressing  and  expanding  pistons,  respectively. 

In  this  phase,  a  shock,  similar  in  strength  to  a  one-dimensional 
shock,  will  form  in  front  of  A  and  surround  the  front  of  the 
cylinder.  Its  strength  decreases  as  the  9=n/2  lire  is  approached; 
then  the  shock  degenerates  into  a  characteristic  surface  surrounding 
an  expansion  zone  in  the  back  of  the  cylinder. 

At  the  same  time,  however,  a  compression  wave  is  sent  by  A 
all  along  the  cylinder  towards  the  rear  and  an  expansion  wave  is 
sent  by  B  towards  the  front.  As  a  result,  the  pressure  decreases 
on  the  front  of  the  body  and  increases  on  its  back.  The  pressure 
distribution  around  the  cylinder  in  the  first  phase  of  the  motion 
is  qualitatively  shown  in  Fig.  2a.  Soon  after  the  motion  starts, 
the  compressed  and  expanded  regions  are  practically  antisymmetrical 
with  respect  to  the  8=rr/2  line.  The  shadowed  region  corresponds 
to  CpX).  The  pressure  waves  did  not  coalesce  yet  into  a  shock. 
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Fig.  2 

The  almost  antisymmetric  pattern  evolves  in  time  along  the 

HAEBI  line  as  shown  in  Fig.  2b. 

The  curvature  of  the  walls  acts  against  the  effects  described 
above.  As  soon  as  the  acceleration  vanishes,  that  is,  at  t*ti, 
the  pressure  at  A  starts  dropping  and  the  pressure  at  B  starts 
increasing.  The  effect  is  shown  in  Fig.  3  by  the  closing  of  some 
of  the  isobars.  Note  that  the  perturbations  produced  at  J  and  L 
(that  is,  the  points  A  and  B  of  the  circle  at  time  tj )  propagate 
along  the  lines  JK  and  LP,  almost  parallel  to  AM  and  BN, 
respectively.  The  effect  above  is  still  antisymmetric.  However, 
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another  affect  takaa  plica  meanwhile,  The  compressions  generated 
along  AJ  and  the  expansions  generated  along  BL  propagate  along 
the  cylinder.  As  a  consequence,  the  pressures  to  the  right  of 
the  line  AQ  tend  to  increase,  and  the  pressures  to  the  right  of 
the  line  BR  tend  to  decrease.  The  isobars  tend  to  bulge,  as 
shewn  in  Fig.  4.  The  effect  is  different  in  the  front  and  bach 
of  the  cylinder  since  the  perturbations  propagate  at  different 
speeds  in  opposite  directions.  The  compression  is  moving  faster 
than  the  expansion;  thus,  the  changes  in  the  isobar  pattern  is 
more  impressive  in  the  back  side  of  the  cylinder.  The  expanded 
zone  tends  to  move  forward,  twisting  the  Cp=0  line  out  of 
symmetry,  and  a  new  compressed  zone  develops  near  point  B  (the 
shadowed  region  in  Figs.  5  and  6  ). 

In  Fig.  6,  note  that  most  of  the  outer  boundary  is  unaffected; 
indeed,  it  moves  outwards  at  the  speed  of  sound  and  cannot  be 
reached  by  perturbations  generated  later  on,  other  than  fast- 
moving  shocks. 

In  a  steady  state,  at  VQ«  1,  (practically  incompressible 

flow) 

(13)  Cp=  A  (2  cos2  0-  — —y  ) 

with  r,  0  polar  coordinates  with  respect  to  the  center  of  the 
cylinder  and  the  x-axis.  The  isobars  appear  in  a  form,  symmetrical 
to  the  0=tt/2  line,  shown  in  Fig.  7.  The  discussion  above,  and 
particularly  Figs.  5  and  6,  show  how  the  transition  from  the 
antisymmetric  pattern  of  Fig.  2  to  the  symmetric  pattern  of  Fig.  7 

can  occur,  in  the  vicinity  of  the  body. 
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MAX.  COMPRESSION 


Zt  rouat  b«  claarly  undaratood,  howavar,  that  the  conplata 
pattarn  of  Fig.  7  (that  is,  axtanding  to  infinity)  ia  thaoraticallv 
achieved  only  after  an  infinita  tirna.  Indeed,  the  pattern  of 
Fig.  7  in  the  vicinity  jf  the  body  ia  aurrounded  by  an  outward 
moving  region  which  containa  the  remains  of  the  original  compression 
and  expansion.  As  pointed  out  above,  no  further  perturbation  can 
reach  and  modify  the  initial  wave  front.  Thua,  at  a  finite  time 
t  a  pattern  similar  to  the  one  shown  in  Fig.  8  should  appear. 

Typical  times  and  distances  can  be  evaluated  approximately 
by  elementary  means.  Since  the  shock  is  weak,  its  velocity  is 
practically  the  same  as  for  a  sound  wave.  Therefore,  the  outer 
boundary  of  the  perturbed  region  can  be  found  as  the  envelope  of 
the  perturbations  issued  by  the  body  points  at  t=0.  The  envelope 
is  defined  parametrically  by 

(14)  [x+xq  (t)-cos(j)]8  +  fy-sincp]8  =  Yts 

if  one  assumes  that  the  speed  of  sound  in  the  perturbed  region 
equals  the  speed  of  sound  of  the  gas  at  rest  and  neglects  the  flow 
velocity.  Elimination  of  cp  between  (14)  and  its  derivative  with 
respect  to  cp  yields  the  »quation  of  the  envelope: 

(15)  [x+xo(t)  ]3  +  y?  =  (l+yTt)8 

which  is  a  circle  centered  at  x=  -xQ(t),  y=0  and  with  a  radius 
equal  to  1  + 

In  a  similar  way,  and  under  the  same  assumptions,  the  envelope 
of  the  perturbations  issued  by  the  body  points  at  time  t*  (when 
the  body  starts  moving  at  a  constant  speed)  is  found  to  be 

(16)  [x+xo(t)-xo(ti)  ]s+ya=  [IVY*  (t-tj)]8 
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that  is,  a  circle  centered  at  x»  -xQ  (t) +xq  ( tx ) «VQ  ( t-tx ) ,  y*0  and 
with  a  radius  equal  to  1*)./y  (t-t:).  The  two  circles  defined  by 
(15)  and  (16)  are  shown  in  Fig.  8:  they  are  marked  FCD  and  GHL, 
respectively.  The  heavy- shadowed  regions  at  the  front  and  rear 
of  the  perturbed  flow  may  contain  pressure  oscillations  produced 
by  compressive  and  expansive  waves  of  decreasing  strength, 
traveling  almost  periodically  along  the  body  as  a  result  of  the 
first  compression  and  expansion;  every  time  one  of  such  waves 
arrives  at  either  A  or  B,  it  sends  a  signal  towards  F  or  D, 
respectively,  and  these  signals  are  alternate  compressions  and 
expansions  of  decreasing  strength. 

Let  us  consider  now  the  case  in  which  the  cyii^^^r  accelerates 
to  a  higher  speed, so  high  that  a  supersonic  bubble  is  expected 
to  appear  at  the  top  of  the  cylinder  in  the  steady  state.  It 
is  well-known  experimentally  that  the  return  from  supersonic  to 
subsonic  flow  tends  to  take  place  abruptly,  across  a  shock  wave. 

The  symmetry  of  the  flow  cbout  the  0=n/2  line,  as  in  Fig.  7,  is 
thus  destroyed  when  a  supersonic  region  exists.  We  can  see  how 
this  happens  by  analyzing  again  the  evolution  of  the  flow  in  time. 

It  is  evident  from  *7igs.  5,6,  and  8  that  the  transition  from 
the  initial,  anti symmetrical  pattern  to  the  steady  state  pattern 
is  characterised  by  a  cor  pression  region  growing  about  point  B 
and  moving  towards  the  front  of  the  cylinder.  If  high  pressure 
differences  exist  between  E  and  B,  the  forward  moving  pressure 
wave  may  coalesce  into  a  shock,  which  will  continue  moving  towards 
E  until  it  reaches  a  position  of  equilibrium.  In  this  case,  Fio.8 
is  replaced  by  Fig.  9,  where  the  imbedded  shock  is  shown,  starting 
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SECTION  III 


OUTLINE  OF  A  NUMERICAL  TECHNIQUE 

The  numerical  computation  is  performed  in  a  region  bounded 
by  the  cylinder  and  by  an  external  boundary  which  is  (except  in 
the  first  phase  of  the  motion)  the  line  separating  the  perturbed 
flow  by  the  outer  flow  at  rest,  with  respect  to  tne  fixed  frame. 
The  numerical  technique  is  composed  of  the  following  parts: 

1.  Computation  of  interior  points. 

2.  Computation  of  points  on  the  cy.linder, 

3.  Computation  of  points  on  the  outer  boundary. 

4.  Test  for  coalescence  of  pressure  waves. 

5.  Imbedded  shock  fitting. 

6.  Treatment  of  the  initial  phase  of  the  motion. 


SECTIOU  IV 
INTERIOR  POINTS 


To  compute  the  interior  points,  the  equations  of  motion  for 
an  inviscid,  compressible  flow  are  written  in  a  polar  frame,  whose 
origin  is  located  at  the  center  of  the  cylinder  and  whose  x-axis 
coincides  with  the  x-axis  of  the  fixed  frame.  Thus,  the  equations 
of  motion  are 


(17) 


P.+uP  +  —  P.+  vu  +  V-ty  —  =  0 
t  r  r  0  r  r  0  T  r 

u  +uu  +  7  u  +  7P  -  —■  +  V  cos  0=0 
w  r  r  w  r  r  o 

VUV  7  v  7  V  ¥  -  vins'° 


St+uSr+  7  V° 
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Here,  u  and  v  are  the  radial  and  transverse  velocity  components 
in  the  moving  polar  frame  (that  is,  relative  to  the  circle), 

P»ln  p,  7»p/p,  S»P-y  In  p. 

In  order  to  maintain  a  certain  fineness  of  the  mesh  near  the 
cylinder  when  the  outer  boundary  moves  away,  a  stretching  of  the 
r-coordinate  is  needed.  The  stretching  is  achieved  by  first 
normalizing  the  r-coordinate  by  means  of  the  transformation 


(18) 


c- 


r-1 

c-1 


where  r=c(q,t)  is  the  polar  equation  of  the  moving  boundary.  Then 
a  new  set  of  variables  (X,Y,T)  is  defined  by 


\  T=  t 

where  a  is  a  function  of  time.  Note  that,  at  the  body,  £=0  and 
X=0,  and  at  the  outer  boundary,  r=l  and  X=l.  If  the  values  of 
X  are  all  equally  spaced,  with  increasing  c  the  values  of  (  tend 
to  accumulate  in  the  vicinity  of  £-0  and  £=1. 

For  any  function  f. 


(20>  W,<r*V-WcV  W'x'Yt+Vt' 

In  turn. 


(21) 


2  tanh 


*r  c-1  '  0“  _^rC9#  t=  “^rct 

(a/2)  _ 1 _  x  1-2X.2>1  l-tanha  ( c/2) 

I- (2  c-1)  s  tanh8  (c/2)  '  c"  2o  2p  1- ( 2'-l j'tanti*  ( a/2) 
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Th«  aquations  of  motion  in  ths  <X,Y,T)  spies  thus  bsconis 
Pt+APX+BPY+Y  (Cu^-Dvy+Ev^+F)  «0 

ut+Aux+b  uy"k3  VH“  0 


(22) 


where 


(23) 


vt+Avx+Bvy+jpx+lpy+0“° 


St+AS  +BSy=0 


c-  xrcr.  D-  -  -  re' ferry  •  E-?xcCn 
F=  7  -  6-  ?C,  J-  7  Sxcc9,  L=  -  Ay 

UV 


B=  -  — ,  H= - —  +  xq  cos  q ,  Q=  — - xQsin  9 

A-  X£Cr[-Cct+u-  |  Cc9]  +Xaat 

To  solve  (22)  numerically,  a  second-order  accurate  scheme 

is  used.  In  the  present  case,  I  abandoned  the  scheme  used  in 

all  my  previous  work  in  favor  of  a  scheme  recently  used  by 
8 

MacCormack  .  Numerical  experiments  performed  on  one-dimensional 
problems  and  on  the  blunt-body  problem  show  that  MacCormack' s 
scheme  yields  the  same  accuracy  with  a  simpler  coding.  The 
advantage  is  particularly  strong  in  vincous  flow  problems  and 
there  are  some  minor  advantages  when  imbedded  shocks  are  present. 
There  are  no  sizeable  reductions  in  computational  time,  however, 
since  the  simpler  scheme  requires  a  double  computation  at  each 
time  step  and  this  is  almost  equivalent  to  computing  once  the 
double  set  of  first  and  second  derivatives. 

In  the  present  scheme,  a  first  set  of  intermediate  values 

—  ^+1 

f  n  m  at  t=(k+l)At  is  computed  by  the  general  formula 
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I 


(24) 


f*+la 

n,m 


fk  + 
n,m 


(A 


rn+r 


-f' 


m 


AX 


n,m  + 


B 


fk  -fk 
n.m+1  n.m 


AY 


+  c)at 


where 


(25) 


•p' 

‘A 

YC 

yE 

O' 

u 

.  A-  - 

G 

A 

0 

0 

V 

9  • " 

J 

0 

A 

0 

.s. 

Lo 

0 

0 

A. 

“B 

0 

-yd 

O' 

rvF1 

0 

B 

0 

0 

H 

L 

0 

B 

0 

9  C®  " 

Q 

0 

0 

0 

0 

0 

and  the  elements  of  A,B  and  C  are  computed  at  the  point  X=nAX, 

k+1 

Y=m/\Y.  The  final  values  f  are  computed  by  the  formula 

n  $  n\ 

fk  J -f-  k+i  7k+1_  fk+l  ?k+l  Yk+1 

{26)  fk+1  =  x-  2*2-  +  (A  +  p-  n^n“  +  C)4? 

n,m  2  AX  AY  2 

where  the  elements  of  A,  B  and  5  are  based  on  the  values  f  k+1  . 

n,m 

It  may  be  noted  that  the  equations  of  motion  are  not 

formulated  in  conservation  form.  Let  it  be  clearly  stated  that 

I  never  used  the  conservation  form  in  any  of  my  previous  works 

9  10 

either,  since  I  do  not  see  any  strong  reason  for  it  ' 


SECTION  V 
BODY  POINTS 

The  body  geometry  in  the  present  problem  is  so  simple  chat 
the  body  points  can  be  computed  by  using  an  integration  scheme 
similar  to  the  one  used  for  interior  points.  Obvious  simplifications 
follow  from  the  vanishing  of  u,X  and  X  .  Only  one  X-derivative, 

Cl 

uv  is  left  in  the  equations.  It  is  advisable  to  approximate  it 
numerically  by  a  three-point  formula  in  both  stages  of  the 
integration  scheme.  Better  results  seem  to  be  obtained,  however, 
if  the  method  outlined  in  Ref.  10  is  used.  A  Cartesian  frame  of 
reference  is  used  at  each  body  point?  its  §-axis  is  normal  to 


1 


IS 


the  wall  and  its  r-axis  is  tangent  to  the  wall.  With  respect  to 
the  cylinder,  the  (?,r|)  frame  has  a  translational  motion  defined 
by  a  constant  velocity,  equal  to  the  flow  velocity  at  the  body 
point,  at  time  (k+l)At.  Let  u  and  v  be  the  components  of  the  flow 
velocity,  relative  to  the  cylinder,  in  the  frame.  A 

characteristic  equation  in  the  ( ^ , t )  plane  (which  moves  along  with 
the  "-axis,  normally  to  itself)  is  written: 

(27)  0=  u .  -  —  (P-P#)-[av  +  x*  cos  91  At 

w  Y  w  p  O 

where  the  values  denoted  by  *  are  computed  at  t=kAt  and  at  a  point 
A*  defined  by 

(28)  (u-a)At,  n*=  vAt 

and  a,  in  (27)  and  (28),  is  the  average  speed  of  sound  between 
the  point  to  be  computed  and  A*.  Note  that  the  zero  in  the 
left-hand  side  of  (27)  represents  the  vanishing  u  component  at 
the  body  point. 

The  values  at  A*  must  be  interpolated  between  body  points 
and  interior  mesh  points  at  time  t=kAt.  A  linear  interpolation 
seems  to  be  sufficient. 


SECTION  VI 

COMPUTATION  AT  THE  OUTER  BOUNDARY 

For  practical  reasons  the  numerical  analysis  is  started  at 

t=0  over  a  region  bounded  by  the  rigid  circle  and  another  circle, 

centered  at  x»  -x  (t  )  and  having  a  radius  equal  to  ‘v— t  ,  where 

o  o  o 

t  «  1.  Such  a  circle  is  the  wave  front  of  the  perturbation 
initiated  at  t=0,  which  propagates  in  all  directions  with  the 
speed  of  sound  of  the  gas  at  rest,  *Jv~,  and  is  left  behind  by  the 
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moving  circle.  So  long  «•  t<tQ,  the  outer  boundary  is  considered 
fixed  and  the  values  of  the  physical  parameters  on  it  are  the 
valu  ,s  of  the  gas  at  rest  (the  velocity  components,  obviously, 
are  not  zero  since  they  are  relative  to  the  moving  circle) . 

After  t=tQ,  the  perturbed  region  must  be  allowed  to  spread 
outside  of  the  initial  boundary.  From  that  moment  on,  the  outer 
boundary  is  then  considered  as  a  moving  shock. 

Obviously,  in  its  rear  part  such  a  line  is  not  a  shock  at 
all  but  a  sound  wave,  which  travels  with  a  velocity  -  xQC08cp 
in  the  normal  direction  at  each  point  (here  cp  is  the  angle  between 
the  normal  to  the  wave  and  the  x-axis).  In  the  first  phase  of 
the  motion  mentioned  above,  the  rear  part  of  the  "shock”  actually 
is  the  leading  edge  of  an  expansion  wave.  However,  the  values 
of  the  physical  parameters  on  the  outer  boundary  and  the  velocity  of 
the  latter  can  be  computed  by  considering  the  boundary  as  a  shock 
all  around  since  a  sound  wave  is  an  infinitely  weak  shock  and 
even  a  weak  expansion  can  be  considered  as  a  shock  of  negative 
strength,  to  within  an  accuracy  of  the  third  order  in  M*-l 
(where  Mr  is  the  normal  component  of  the  Mach  number,  relative 
to  the  shock ) . 

The  computation  of  the  shock  proceeds  as  outlined  in  Ref.  11. 
The  three  Rankine-Hugoniot  conditions,  the  conservation  of 
tangential  velocity,  and  one  characteristic  equation  written  for 
the  interior  of  the  computational  region  are  sufficient  to 
determine  the  shock  velocity  and  the  values  of  P,u,v,  S  behind 
the  shock.  A  Cartesian  frame  of  reference  is  used  at  each  shock 
point;  its  *-axis  is  normal  to  the  shock  and  its  ir-axis  is  tangent 
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to  ths  shock,  as  it  is  st  tins  t-  (k+l)At.  With  rsspsct  to  ths 
cylinder,  ths  ({,n)  frsns  has  a  translational  notion  defined  by 
a  constant  velocity,  equal  to  the  tangential  canponent  of  the 
flow  velocity  at  the  shock  at  tine  (k+l)At.  Let  u  and  v  be 
the  components  of  the  flow  velocity,  relative  to  the  cylinder, 
in  the  (?,r0  frame.  Let  N1#  Na  and  T1,Ta  be  the  directive  cosines 
of  the  5-axis  end  the  r^-axis,  respectively,  with  respect  to  the 
r-axis  and  the  normal  to  it,  as  they  are  at  t"(k+l)At  for  the 
shock  point  to  be  computed.  Then  the  tangential  component  of 
the  velocity  at  the  shock  point  is 

(29)  vs=  -xo(T1cos0-rasin0) 

The  normal  component  of  the  velocity  at  the  shock  point, 

in  front  of  the  shock  is 

(30)  u„  =  -x  (Mi cos  0-Nasino) 

S  \  o 

The  normal  component  of  the  velocity  relative  to  the  shock, in 
front  of  it,  is 


(31) 


V  .«  u  -W 
n  rel  Si 


where  W  is  the  shock  velocity  in  the  normal  direction*  With 
Vnrel  a  relative  normal  Mach  number  is  made. 


(32) 


M*  Vn  rel  /S 


and  the  first  Rankine-Hugoniot  equation  yields  the  relative  normal 

e  shock, 

1+  M* 


component  of  the  velocity  behind  the  shock, 

t- 1 


(33) 


n  rel  n  rel^  y+1  M* 
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The  second  and  third  Rankine-Hugoniot  aquationa  yiald  tha  danaity 
and  praaaura  bahind  tha  shock, 


(34) 

.  n 

n  rel 

(35) 

P_  Iv-l)n-(y-l) 

V+1-  ( Y-l)  p 

From  these  equations,  the  values  of  P  and  S  behind  the  shock  are 
found: 


(36)  P=  In  p 

(37)  S-  P-yln  p 


A  characteristic  equation  in  the  (?,t)  plane  (which  moves  along 
with  the  g-axis,  normally  to  itself)  is  written: 

(38)  u  =u*+  ^(P-P*)+ [a*v  -V  (NiCOse-»3sine) lAt 

S  Y  r*  ° 

where  the  values  denoted  by  *  are  computed  at  t=kAt  and  at  a 
point  A*  defined  by 


(39)  |ua-a  |  At  ,  rf=  -vg  At 

and  a,  in  (38)  and  (39),  is  the  average  speed  of  sound  between 
the  point  to  be  computed  and  A* .  Since,  on  the  other  hand. 


(40) 


u  -V  .+  W 
s  n  rel 


the  set  of  equations  above  can  be  used  iteratively  to  determine 
all  the  unknowns  contained  therein. 

The  values  at  A*  must  be  interpolated  between  shock  points 
•  and  interior  mesh  points  at  time  t«kAt.  The  interpolation  can 
be  coded  in  many  different  ways.  A  well- formed  shock,  followed 
by  a  rather  flat  distribution  of  values,  presents  no  difficulties. 
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However,  for  a  shock  generated  by  coalescence  of  compression  waves, 

the  distribution  of  values  behind  it  is  very  steep.  A  simple 

linear  interpolation  tends  to  underestimate  the  pressure  at  A*. 

Consequently,  the  computed  value  of  P  at  the  shock  is  too  low 

and  the  shock  does  not  build  up  strength.  The  errors  generated 

at  the  shock  propagate  inside  the  computational  region,  producing 

oscillations  which  are  not  physically  justified.  This  aberration 

9 

has  been  detected  in  one-dimensional  problems  and  I  suggested 
a  different  interpolation  scheme  which  should  work  well  under 
any  circumstances.  It  should  be  noted  that,  if  the  linear  scheme 
is  used  in  a  one-dimensional  problem,  the  wiggles  tend  to 
disappear  once  the  piston  reaches  a  constant  speed  (because  the 
piston  keeps  sending  high-pressure  signals  towards  the  shock, 
thus  giving  it  a  chance  to  gain  strength,  eventually  ) .  Therefore, 
the  aberration  is  disturbing  but  not  lethal.  In  the  present 
problem,  instead,  the  pressure  behind  the  shock  must  decrease, 
and  strongly,  in  the  second  phase  of  the  motion  described  above. 

The  shock  cannot  correct  its  speed  and  shape,  and  the  errors 
generated  at  the  shock  tend  to  grow  bigger,  inside  the  computational 
region,  as  the  overall  pressure  decreases.  The  oscillations 
may  eventually  grow  beyond  control.  It  is  thus  imperative  to 
use  a  better  interpolation  scheme  in  the  shock  computation; 
therefore,  the  scheme  suggested  in  Ref.  9  has  been  adapted  to 
the  present  two-dimensional  problem. 
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SECTION  VII 


RADIAL  STRETCHING 

To  complete  the  outline  of  the  numerical  scheme,  a  word  must 
be  spent  on  the  definition  of  a  in  (19)  .  If  c«  1,  the  relationship 
between  a.  and  X  is  almost  linear.  Therefore,  equal  intervals 
along  the  X-axis  correspond  to  almost  equal  intervals  between 
shock  and  body.  With  increasing  distance  between  shock  and  body, 
the  mesh  size  in  the  physical  plane  tends  to  grow  too  large. 

Now,  there  are  two  regions  in  the  flow  field  where  we  must 
maintain  a  high  radial  resolution.  One  is  a  ring  surrounding 
the  body,  where  we  want  to  achieve  high  accuracy  since  the  values 
at  the  body  are  the  most  important  from  a  practical  standpoint. 

The  other  is  the  ring  (FCDGHL)  shown  in  Fig.  8,  where  the 
pressure  changes  rapidly  from  low  to  high  values  and  to  low  values 
again  for  0>  n/2  and  has  a  similar,  but  reversed,  behavior  for 
9<  tt/2  .  Such  changes  are  due  to  the  fact  that  the  initial 
perturbations  (a  strong  compression  coalescing  into  a  shock  and 
a  strong  expansion,  respectively  for  0>  tt/2  and  9<  n/2)  are 
followed  by  perturbations  of  the  opposite  kind.  The  thickness 
of  the  ring  is  of  the  order  of  the  time  spent  to  accelerate  the 
body,  multiplied  by  the  speed  of  sound,  that  is,  Jy  tt/uj.  Therefore, 
the  nodes  must  concentrate  in  the  vicinity  of  X^O  and  X«1  and 
this  justifies  the  choice  of  the  stretching  function  (19) .  In 
addition,  a  should  start  increasing  as  soon  as  the  first  row  of 
interior  points  reaches  a  distance  d  from  the  body  which  is 
considered  as  the  greatest  to  be  allowed  for.  From  (18)  and  (19), 
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and  forcing  r-l-d  whan  X-AX,  wa  obtain 

(41)  .  ^[1+  tanha(AX-.5)  j 

tanh(o/2) 

Than 

(42) 

d  2  tanh (a/2) 

at*  “  (c-l)B  ct  (ax- .5) Ll-tanh*  a (AX- .5) J-2  tanh a ( AX- . 5)  cosscha- 

This  equation  car  be  conveniently  used  to  increment  a  at  each 
step  in  order  to  keep  d  practically  constant.  In  (41)  and  (42), 
the  values  of  c  and  c^.  are  assumed  to  be  the  ones  on  a  typical 
radius,  for  example,  on  the  x-axis  in  front  of  point  A. 

SECTION  VIII 

RESULTS  FOR  A  FULLY  SUBSONIC  FLOW 

To  better  understand  how  t.»e  flow  develops  in  a  simple 
case,  where  no  imbedded  shocks  are  present,  I  analyzed  the  case 
in  which  the  Mach  number  of  the  uniform  flow  at  infinity, 
relative  to  the  circle,  is  equal  to  0.1,  and  the  circle  accel¬ 
erates  slowly  («>=tt),  making  an  attempt  to  perform  an  exhaustive 
study.  A  preliminary  computation  was  made,  with  a  simplified 
stretching  in  the  radial  direction.  The  value  of  d  was  taken 

as  0.05,  and  t  =  0.2.  The  mesh  had  7  intervals  in  the  radial 
o 

direction  and  18  intervals  circumferentially.  Some  of  the  resulting 
isobars  are  presented  in  Fig.  10.  They  show  qualitative  agreement 
with  the  pattern  foreseen  in  Fig.  5.  However,  early  wiggles  indicate 
a  poor  resolution  of  the  mesh.  From  a  more  detailed  plot,  which 
cannot  be  reproduced  in  a  small  scale,  it  appears  that  more  and 
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more  wiggles  tend  to  develop  it  the  beck  of  the  cylinder  towards 

* 

the  last  computed  step  (at  t»5.25,  after  300  cycles)  . 

Increasing  the  number  of  radial  intervale  to  28  and  using 
the  sketchings  defined  above  with  d«0.05  definitely  improves 
the  calculation.  The  isobar  pattern  so  obtained  (Fig.  11)  can 
be  compared  with  Fig.  10.  No  wiggles  appear  and  the  pattern 
follows  closely  the  predictions  of  Section  II.  Note  that  the 
Cp=0  line  issuing  from  the  top  point  on  the  cylinder  almost 
reaches  the  forward  stagnation  point  bringing  a  strong  expansive 

effect,  and  then  recedes  towards  its  steady  state  location,  30° 

*  ★ 

above  the  stagnation  point.  After  650  cycles  ,  the  pressure 
distribution  along  the  circle  is  the  one  to  be  found  in  the 
steady  state  to  within  1%.  Figs.  12,  13  and  14  show  in  a  different 
way  how  the  steady  state  about  the  circle  is  reached.  Fig.  12 
is  a  plot  of  isobars  in  the  plane  of  motion  at  the  400th  cycle 
( t=5 . 06) .  This  is  about  the  time  when  the  forward  stagnation 
point  undergoes  the  maximum  expansion.  Fig.  13  is  a  similar  plot 
at  the  700th  cycle  (t*10.36).  The  steady  state  is  achieved 
within  a  circle  centered  at  the  center  of  the  body  and  having  a 
radius  twice  the  radius  of  the  body. 

Finally,  Fig.  14  is  a  similar  plot  at  the  1000th  cycle 
(t=15.66)+.  The  region  within  a  circle  centered  at  the  center  of 
the  body  is  practically  in  a  steady  state.  The  result  is 
consistent  with  the  propagation  of  signals  at  a  speed  practically 
equal  to  the  speed  of  sound  in  the  gas  at  rest:  if  the  steady 
* 

One  minute  of  CDC  6600  time. 

** 

About  6  minutes  of  CDC  6600  time. 

+  About  eight  minutes  of  CDC  6600  time. 
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Figure  12 
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Figure  13 
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Figure  14 

state  front  is  at  r=2  when  t=10  and  travels  at  a  speed  equal  to 
Jy  =  1.183,  it  reaches  r*8.7  at  t*15.6. 

The  history  of  the  pressure  distribution  on  the  symmetry  line 
in  front  of  the  body  is  given  by  Fig.  15.  The  pressure 
coefficient  is  plotted  as  a  function  of  r  (the  body  is  at  the  left 
and  the  initial  perturbation  front  at  the  right)  every  100 
computational  cycles.  One  can  see  how  the  initial,  strong 
compression  tends  to  form  a  shock.  However,  its  tail  is  promptly 
reduced  by  the  following  expansion.  Then  a  moderate  compression 
builds  up  in  the  vicinity  of  the  body  and  the  steady  state  spreads 
outwards.  The  original  compression  and  the  following  expansion 
are  reduced  in  strength  when  the  computation  stops,  but  they  are 
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far  from  having  been  wiped  out. 

Figs.  16,17  and  18  confirm  the  fact  that  a  steady  state  has 
been  reached.  The  exact  values  of  Cp  in  an  incompressible,  steady 
flow  are  given  by 

(43)  C  =  2c°|2  9 - \- 

p  r  r 

In  Figs.  16,  17  and  18  the  solid  lines  represent  the  exact  values 
(43)  on  the  front  and  rear  stagnation  lines  and  on  the  o*tt/2  line. 
The  computed  values  at  the  700th  cycle  are  plotted  in  Fig.  16  and 
the  computed  values  at  the  1000th  cycle  are  plotted  in  Fig.  17. 

Not  only  is  it  clear  thac  the  steady  state  grows  in  size  as 
explained  above,  but  the  accuracy  of  the  results  is  proved.  In 
a  similar  way,  the  solid  line  in  Fig.  18  represents  the  exact 
values  (43)  along  the  body  (r=l)  and  the  computed  values  are 
plotted  for  a  comparison. 

Fig.  19  shows  the  pressure  on  the  front  stagnation  line 
immediately  behind  the  perturbation  front,  as  a  function  of  time. 
It  appears  that  the  shock  starts  building  up  at  the  time  predicted 
by  the  one  dimensional  theory  if  the  piston  moves  with  the  same 
law  as  the  present  body. 

Let  is  be  noted  explicitly  that  the  stretching  of  the  radial 
coordinate  and  the  number  of  nodes  are  necessary  and  strictly 
sufficient  to  get  the  results  presented  above,  which  I  consider 
satisfactory  from  both  a  theoretical  and  a  practical  standpoint. 

If  fewer  nodes  are  taken  in  the  outer  ring,  wiggles  appear  in  the 
far  field  almost  simultaneously  with  the  beginning  of  the  steady 
flow  in  the  near  field.  The  wiggles  tend  to  move  inwards  as  time 
increases  and  eventually  they  perturb  the  entire  steady  state 
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Figure  18 


region. 

In  conclusion,  we  can  issue  a  statement  about  the  minimum 
computational  time  needed  for  determining  the  steady  state  at  the 
surface  of  the  body.  Starting  at  the  instant  at  which  the  ac¬ 
celeration  of  the  cylinder  stops,  one  should  let  the  expansion 
wave  travel  all  the  way  from  the  rear  stagnation  point  to  the 
fore  stagnation  point  and  then  again  travel  all  the  way  back  to 
the  rear  stagnation  point.  The  slowest  speed  of  the  wave  is 
given  by  the  difference  between  the  speed  of  sound  and  the  maximum 
modulus  of  the  tangential  velocity  component.  The  latter  is 
reached  at  9=  90°  and  is,  according  to  the  Rayleigh- Jantzen 
second  order  approximation  for  compressible  flow  (Ref.  14  ) 
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Figure  19 

The  steady  state,  thus,  starts  building  up  at  the  surface 
of  the  cylinder  approximately  at 

(44)  t=  1 1  +  - —7 - 

-g  M®8 )  ] 

if  the  speed  of  sound  around  the  cylinder  is  considered  equal  to 
the  speed  of  sound  at  infinity. 

In  the  present  case,  M  =  .1,  ti  =  .5  and  tz7.2.  With  a  safety 

00 

factor  of  1.5,  one  can  safely  interrupt  the  computation  at  t=10, 
that  is,  after  700  computational  cycles  (about  6  minutes  of  CDC 
6600  time) . 

For  a  cylinder  one  meter  in  diameter  (*ref~  *5)  accelerating 
to  M*  .1  in  a  standard  atmosphere  at  sea  level  (8^360  m/sec), 
tref=  .00164.  Then  the  physical  time  corresponding  to  t=10  equals 
.0164  sec.  The  ratio  t  between  computational  time  and  real  time 
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is 


.  360 

C"  .0164 


(45) 


-  22,000 


In  the  present  computations,  the  safety  coefficient  applied  to 
the  Courant-Friedrichs-Lewy  stability  formula  was  taken  substantially 
lower  than  necessary,  and  the  size  of  the  initial  region 
substantially  smaller  than  the  distance  of  formation  of  the  shock 
wave.  By  relaxing  both  conditions,  one  should  be  able  to  reduce 
the  computatior  time  to  about  4  minutes,  thus  reducing  t  to  about 
15,000.  This  .eems  to  be  a  lower  boundary  for  t,  with  the  present 
generation  of  computers. 
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SECTION  IX 


A  CRITICAL  SURVEY  OF  OTHER  APPROACHES  TO  THE  PROBLEM 

At  this  stage,  we  have  a  complete  analysis  of  the  transition 
from  the  incipient  motion  to  a  steady  state  regime  and  a 
consistent  matching  of  the  physical  evolution  with  its  numerical 
description.  Such  a  physically  oriented  study  may  also  help  us 
in  judging  the  reliability  of  other  numerical  treatments  of  the 
steady  state  problem,  all  generally  ascribed  to  the  category  of 
"time-dependent  calculations". 

A)  The  flow  relative  to  the  cylinder  is  studied  in  a  frame 
fixed  to  it.  An  external  boundary,  also  fixed  with  respect  to 
the  cylinder,  is  chosen  arbitrarily.  It  is  assumed  that  the  flow 
parameters  on  such  a  boundary  are  those  of  the  uniform  flow  at 
infinity.  An  arbitrary  initial  distribution  of  values  is  assumed 
in  the  computational  region  and  the  unsteady  flow  is  computed. 

Such  an  approach,  which  seems  to  have  been  popular  some  years  ago, 
was  supported  by  the  belief  that,  since  the  external  boundary  is 
taken  at  a  respectable  distance  from  the  body,  the  conditions  on 
it  are,  to  all  practical  purposes,  the  same  as  at  infinity  (see, 
for  example.  Ref.  5) .  I  have  expressed  my  disagreement  in  Ref. 

12.  Here  I  wish  to  strengthen  it  on  the  light  of  the  analysis 
above.  We  have  seen  that  the  steady  state  results  from  the 
interplay  of  signals  travelling  around  the  body  and  outwards .  We 
have  seen  the  disastrous  effects  on  the  near  field  of  perturbations 
proceeding  inwards  from  the  far  field.  Now,  in  the  first  phases 
of  a  computation  like  the  one  described  in  this  section,  the  signals 
from  the  body  travel  towards  the  outer  boundary  and  cannot  get  out. 
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The  boundary  sends  signals  inwards,  which  have  very  little  to  do 
with  the  outgoing  waves.  Indeed,  such  a  boundary  does  not  act  as 
a  rigid  wall  but  it  does  not  act  as  an  infinite  capacity  where 
the  pressure  is  constant  since  velocity  and  entropy  are  also 
kept  constant  on  it.  Consequently,  one  may  expect  to  find  an 
initial  stage  in  which  the  results  in  the  near  field  tend  to 
approach  the  steady  state  solution,  followed  by  a  second  stage  of 
increasing  inaccuracy. 

A  computation  of  this  kind  (assuming  local  consists  cy  of  the 
numerical  scheme  with  the  equations  of  inviscid  unsteady  motion) 
is  surely  a  time-dependent  one.  However,  it  does  not  converge  to 
the  desired  steady  state  pattern  since  it  describes  a  problem  with 
boundary  conditions  different  from  the  ones  consistent  with  the 
time-dependent  evolution  of  a  flow  extending  to  infinity.  In 
certain  cases,  the  results  are  deceptively  smooth.  This  is  due 
to  the  presence  of  artificial  viscosity  in  the  numerical  scheme. 
Clearly,  there  is  no  reason  for  the  results  of  a  dissipative 
calculation  to  be  consistent  with  the  time  solution  of  an  inviscid 
problem. 

B)  The  infinite  field  about  the  obstacle  is  mapped  onto  a 
finite  computational  region  by  means  of  a  stretching  of  coordinates 
in  the  radial  direction.  Again,  arbitrary  initial  conditions  are 
assumed.  I  proposed  this  technique  in  Ref.  12  and  some  applications 
were  made  in  Ref.  7.  By  so  doing,  uniform  values  can  be  imposed 
on  the  outer  boundary  since  the  latter  is  actually  at  an  infinite 
distance  in  the  physical  plane,  and  in  principle,  it  is  reached 
by  outgoing  waves  only  at  t=®.  At  the  light  of  the  present  study. 
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it  is  clear  that  a  steady  state  tends  to  build  up  in  the  vicinity 
of  the  body  but  the  computation  eventually  becomes  unstable. 

Indeed,  the  spacing  between  nodes  increases  from  the  body  to 
infinity.  The  resolution  is  extremely  poor  where  it  should  be  very 
high  (within  the  ring  which  follows  the  initial  perturbation 
front).  In  addition,  in  our  exercise  reported  in  Ref.  7,  we 
assumed  impulsive  initial  conditions.  The  initial  phase  of  the 
motion,  thus,  was  characterized  by  very  strong  perturbations. 

The  technique  was  unable  to  handle  shocks.  As  a  consequence,  the 
unsteady  outgoing  ring  was  full  of  numerical  oscillations.  As 
soon  as  these  entered  the  region  of  poor  resolution,  they  started 
feeding  errors  inwards.  The  flow,  which  had  reached  an  almost 
perfect  pressure  distribution  near  the  body  after  1000  computational 
steps,  was  competely  destroyed  at  step  2000.  We  can  see  now  that 
the  instability  of  the  computation  is  explained  again  by 
inconsistency  of  the  numerical  analysis  with  the  physical  pattern. 
The  same  code,  with  initial  conditions  given  by  the  incompressible, 
steady  flow  field,  went  through  2000  cycles  with  no  symptoms  of 
deterioration  and  errors  less  than  0.0002.  It  should  also  be 
noted  that  very  good  results  in  the  vicinity  of  a  body  can  be 
obtained  by  mapping  the  infinite  physical  field  onto  a  finite 
computational  field  and  letting  the  body  st^-t  from  rest,  in  those 
cases  where  no  shocks  build  up  along  the  perturbation  front. 
Obviously,  the  computation  is  to  be  halted  where  the  steady  state 
is  reached  near  the  body  and  before  the  poor  far  field  resolution 
damages  the  near  field.  Problems  of  this  nature  will  be  discussed 
in  a  forthcoming  paper  by  B.  Grossman. 
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C)  Finally,  tha  computation  praaantad  in  Ref.  6  conaiata  of 
aaauming  impulaive  initial  conditiona  and  atarting  at  a  slightly 
later  time.  A  perturbation  front,  in  the  form  of  a  shock  wave, 
similar  to  the  one  described  in  this  report,  is  assumed  around 
the  body  and  let  to  expand.  Unfortunately,  an  initial  flow  field 
must  be  defined  between  shock  and  body.  Such  a  definition  is  bound 
to  be  arbitrary  (in  the  present  study,  inataad,,  the  initial 
conditions  are  not  arbitrary).  The  result^  presented  in  Ref.  6 
are  not  so  detailed  to  allow  an  analysis  as  minute  as  in  this 
report.  However,  in  the  light  of  the  present  analysis,  it  appears 
that  the  initial  abitrariness  is  forever  trapped  within  the 
computational  region  and  it  should  eventually  damage  the  steady 
state  around  the  body  when  the  latter  tends  to  spread  outwards. 

By  proceeding  as  in  Ref.  6,  it  seems  to  be  difficult  to  stop  the 
calculation  at  the  right  time  because  there  is  no  criterion 
available  to  determine  whether  and  when  a  steady  flow  is  reached 
in  the  vicinity  of  the  body.  The  Author  circumvents  the  onset 
of  instabilities  by  using  a  low  order  integration  scheme  and 
numerical  dissipation,  two  devices  which  I  do  not  recommend. 


39 


SECTION  X 


REFERENCES 

1.  Southwell,  R.V.,  Relaxation  Methods  in  Theoretical  Physics, 
Oxford,  1946, 

2.  Hayes,  W.D.  and  Probstein,  R.F.,  Hypersonic  Flow  Theory,  Vol.I, 
page  438,  Academic  Press,  New  York,  1966. 

3.  Crocco,  L.,  A  Suggestion  for  the  Numerical  Solution  of  the 
Steady  Navier-Stokes  Equations,  AIAA  J.,  .3,  1824,  1965. 

4.  Evans,  M.W.  and  Harlow,  F.H.,  Calculation  of  Unsteady  Super¬ 
sonic  Flow  Past  a  Circular  Cylinder,  ARS  J.,  January  1959, 
page  46. 

5.  Magnus,  G.,  Gallaher,  W.,  and  Yoshihara,  H.,  Inviscid  Super¬ 
critical  Airfoil  Theory,  AGARD  Specialist'  Meeting  on 
Transonic  Aerodynamics,  Paris  18-20  September  1968. 

6.  Kentzer,  C.,  Computations  of  Time  Dependent  Flows  on  an  Infinite 
Domain,  AIAA  8th  Aerospace  Sciences  Meeting,  New  York,  1970, 
Preprint . 

7.  Mackenzie,  D.,  and  Moretti,  G.,  Time  Dependent  Calculation 
of  the  Compressible  Flow  About  Airfoils,  AGARD  Specialist' 
Meeting  on  Transonic  Aerodynamics,  Paris,  18-20  September 
1968. 

8.  MacCormack,  R.W.,  The  Effect  of  Viscosity  in  Hypervelocity 
Impact  Cratering,  AIAA  7th  Aerospace  Sciences  Meeting, 

Paper  No.  69-354,  1969. 

9.  Moretti,  G.,A  Critical  Analysis  of  Numerical  Techniques: 

The  Piston-Driven  Inviscid  Flow,  PIBAL  Report  No.  69-25, 

July  1969. 


40 


10.  Moretti,  G . ,  The  Choice  of  Tine-Dependent  Technique  in  Gas 
Dynamics,  PIBAL  Report  No.  69-26,  July  1969. 

11.  Moretti,  G.  and  Abbett,  M.,  A  Tine -Dependent  Computational 
Method  for  Blunt  Body  Flows,  AIAA  J.,  4,  2136,  1966. 

12.  Moretti,  G.,  Importance  of  Boundary  Conditions  in  the  Numerical 
Treatment  of  Hyperbolic  Equations,  Phys.  Fluids,  Supplement 

II,  11-13,  1969. 

13.  Lifshitz,Yu.B. ,  The  Computation  of  Transonic  Flew  Past  a 
Symmetric  Profile,  Mekh.  Dz  i  Gaz.  page  52,  1969. 

14.  Howarth,  L., (editor).  Modern  Developments  in  Fluid  Dynamics, 
High  Speed  Flow,  Page  275,  Oxford  1953. 


41 


Tmrurrri.fTnr 


DOCUMtNT  CONTROL  DATA  -RAD 

f ttturlly  tlatlllltallan  at  III!*,  b*dy  a 1  abtlfrl  and  Induing  annalailan  mu*l  b*  aatarad  »b*n  rh*  *  frail  raparl  1 •  tlaaalllad J 

Polytechnic  Institute  of  Brooklyn 

Dept,  of  Aerospace  Engrg.  4  Appl.  Mechanice 
Route  110,  Farmingdale,  New  York  11735 

Unclassified 

a*.  «mou» 

AIAOAT  TITLI 

TRANSIENT  AND  ASYMPTOTICALLY  STEADY  FLOW  OF  AN 
GAS  PAST  A  CIRCULAR  CYLINDER 

INVISCID,  COMPRESSIBLE 

4.  OKIC  aiativk  noth  (Type  ol  report  on  4  fnc/uefre  defee; 

•  autmoAiIi  (Ft  ft  nemo,  mlddl e  initial,  tart  nemo) 

Gino  Moretti 

•  MC^OMT  OATI  !?•.  TOTAL  NO  OP  »*01t  Tib.  NO  O* 


•  MK^OWT  O  A  T  K  ?•.  TOTAL  NO  O*  PAOIB  Tfc.  NO  O*  AIM 

ADril  1970  41  14 


•«  CON  T  A  AC  T  OA  CAAHT  NO  M.  O  AlO  IN  A  TO  A'l  AIAOAT  NUM  B I  A|ll 

DAHC04-69-C-0077 

b.  MOJIC  T  NO 


PIBAL  REPORT  NO.  70-20 

c  ARPA  Order  No.  1442  |»fc.  OThIA  AIAOAT  NOItl  (Any  othat  number*  that  may  6«  +>md 

thla  report) 

d.  Program  Code  9E30 

10  DIITAIIUTlON  ITATIMINT 

Distribution  of  this  document  is  unlimited. 

11  IUAPlIMCNTAAV  NOTES  =  IJ  IAONIOAIN6  Win  rAAV  ACTIVITY 

U.S.  Army  Research  Of f ice-Durham 
Box  CM 

Duke  Station,  Durham,  N.  Carolina 

O  A  n  t  T  a  A  C  T 

A  numerical  analysis  is  made  of  the  inviscid  flow  produced  by  a 
cylinder  which  accelerates  from  a  state  of  rest  to  a  constant,  subsonic 
speed  in  a  gas  at  rest.  All  features  of  the  numerical  solution  are 
explained  on  physical  grounds.  Consequently,  ways  are  suggested  to 
computed  steady  subsonic  flows  around  obstacles  with  a  maximum  accuracy 
and  a  minimum  computational  time. 


Unclassi fied 

Security  Cl»**ific»ti«n 


DD  /®o?..1473 


Unclassified 


uzsaifflucn 


Numerical  techniques 
Time-dependent  techniques 
Subsonic  flow  past  a  cylinder 


Unclassified 


Security  Clarification 


